% PICCS test
%-----------------------------

%%============================
%get test data ready

% image data load from 1.bmp
img = double(imread('1.bmp'));% 64 * 64 image
img = img(:,:,1);
x = img(:);


n =64;
N = n*n;


%-------phamton test----------------
%img = phantom(n);
%x = img(:);
%-----------------------

% create (sparse) differencing matrices for TV
Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);


%set projection num
 L = 16;
 
%measure matrix
matrix = getSampMatrix(n,L);
A = double(matrix);

% measurements
y = A*x;

%lamda and xp
lamda = 0;
theta = 1:(180/L):180;
[R] = radon(img,theta);
xp = iradon(R,theta,n);
xp = xp(:);


%%reconstruction
cvx_begin
    variable xrec(N)
    %minimize( norm(xrec-xp,1) )  
    minimize(lamda*norm(xrec - xp,1)+ (1-lamda)*norm(xrec,1))
    subject to
        y  == A * xrec ;           
cvx_end
 
figure(2)
subplot(3,1,1);
imshow(img);
title('orig');
subplot(3,1,2);
imshow(reshape(xp,n,n));
title('back projection reconstructed');
subplot(3,1,3);
imshow(reshape(xrec,n,n));
title('reconstructed')
